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Abatract :  Spatially  invariant  point  spread  functions  are  a  common  model 
for  image  and  signal  degradation.  In  general,  the  process  of  reversing 
gaussian  blur  is  unstable,  and  cannot  be  represented  as  a  convolution 
filter  in  the  spatial  domain.  However,  if  we  restrict  the  domain  of 
allowable  functions  to  polynomials  of  degree  no  greater  than  N,  then  an 
inverse  filter  exists.  We  use  Herraite  polynomials  to  represent  kernels 
which  can  be  used  to  deblur  polynomial  data  which  has  been  degraded  by 
a  known  amount  of  gaussian  blur.  For  fixed  N,  the  corresponding  kernel 
gives  stable  deblurrlng  among  the  class  of  functions  which  are  gaussian 
filtered  versions  of  data  well  approximated  by  polynomials  degree  N  and 
less. 
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I.   Introduction 

Realistic  Imaging  and  signal  sampling  systems  necessarily 
Introduce  spatial  degradation  of  the  original  data.  Frequently,  the 
degradation  process  is  modeled  as  a  spatially  invariant  gausslan 
convolution.  Thus  if  f(x)  Is  the  original  data,  x  e  R",  then  the 
observed  data  is 


h(x)  =  i   K(x  -  x',  t)  f (x')  dx'  , 
R" 


where 


K(x.t)   =  1     e-l-1'/^^ 


(ATit)"/2 

is  the  gausslan  kernel,  whose  extent  is  parameterized  by  t  >  0,  and  is 
normalized  to  have  unit  mass.  How  can  the  oiginal  data  £(x)  be 
reconstructed  when  only  h(x)  and  the  amount  of  blurring  t  is  known? 

As  is  well  known,  in  general  the  data  f(x)  cannot  be  reconstructed 
from  h(x).  Not  all  functions  h(x)  arise  as  blurred  versions  of  some 
original  data  f(x),  and  even  if  an  Ideal  function  f(x)  exists, 
arbitrarily  small  inaccuracies  in  the  representation  of  h(x)  can  lead 
to  large  inaccuracies  in  the  reconstruction  of  f(x). 

Nonetheless,  as  is  also  well  known,  deblurring  can  in  practice  be 
accomplished  by  high-emphasis  filtering.  The  inherent  numerical 
instabilities  in  general  cause  no  problems.  This  apparent 
contradiction  is  resolved  by  the  observation  that  the  deblurring 
process  is  stable  when  restricted  to  a  suitable  class  of  functions,  and 
that  this  class  subsumes  most  naturally  occurring  signals. 

In  this  paper,  we  present  kernels  which  can  be  used  to  deblur  a 
fixed   amount   of   gausslan  blur,   and  accomplish  this  inverse  process 
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stably  among  original  functions  f(x)  that  are  well  approximated  by 
polynomials  of  fixed  degree.  These  kernels  will  provide  exact 
reconstructions  of  (sufficiently  low  order)  polynomial  original  data. 
The  existence  and  use  of  these  kernels  has  been  presented  elsewhere 
[John];  our  analysis  is  new  only  in  that  we  emphasize  the  use  of 
Hermite  polynomials,  and  give  explicit  formulas  in  terms  of  Hermite 
polynomials. 

The  problem  of  reconstructing  f(x)  given  h(x)  and  t  is  the  inverse 
heat  equation  problem,  since  the  function  h(x)  represents  a 
distribution  of  heat  after  t  units  of  time,  where  f(x)  is  the  initial 
t  =  0  distribution.  Solving  Che  heat  equation  backwards  in  time  is  an 
interesting  inverse  problem  with  applications  in  image  enhancement, 
signal  and  image  representation,  and  perhaps  even  modeling  of  neural 
processes  for  analyzing  data. 

In  the  next  section,  we  formulate  the  deblurring  problem  in  one 
dimension.  For  data  in  higher  dimensions,  we  shall  simply  appeal  to 
the  separability  of  the  gaussian  kernel.  Our  formulation  is  strictly 
in  a  continuous  Euclidean  domain.  A  separate  analysis  for  discrete 
data  is  possible  (see,  e.g.,  Pratt],  but  presents  different 
difficulties. 
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II.   Formulating  the  Deblurring  Problem 

Consider  the  opertor  Tl^    defined  on  L  (R)  bv  the  equation 


(fi^fXy)  =  ]  -1_  e"'^  /^^  f(y-x)  dx 
-co  Z/tTE 


For   t   >  0,  n^  Is  a  compact  symmetric  bounded  linear  ooerator  on  L  (R) 
mapping  into  L  (R).   This  operator  has  many  special  properties,  such  as 


"t  °  "s  ="t+s 
and 

u(x,t)  =  (n^f)(x)      satisfies   Au  =  u^  ,  u(x,0)  =  f(x); 

see  [Bars,  John,  Schechter] .   If  we  denote  the  Fourier  transform  of   a 
function  g(x)  by  g(a)),  then  fi^.  is  a  multiplier  operator  given  by 

in^ffCoi)   =  e-i'^  ^    f(u)) 

By  means  of  this  formula,  ^      can  be  extended  to  operate  on  the  Schwartz 
class  S^'  of   Fourier   transformable   distributions   [Rudin,   Funotlonal 
Analysirs] «   In  particular,  Q  ^f    is  defined  for  any  polynomial  f. 
We  will  specialize  to  the  case  t  =  1/^,  and  set 


Tf  =  fij/^f  =  -L  e"'^  *  f  .  (2.1) 


By   suitably  scaling   the   spatial  parameter  x  e  R,  f^^  ,  t  >  0,  can  be 
seen  to  be  equivalent  to  T  operating  on  a  rescaled  version  of  f,  i.e., 
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(n^fXy)   =  (Tf)(y/2/t)  , 
where   ?(x)   =  f(2  ft   x). 

Thus  the  Invertlbi llty  of  ^^    is  settled  by  Inverting  T. 
From  the  Fourier  multiplier  formula 

(Tf)'(u))  =  e^  /^  f(a))  ,  (2.2) 

and  the  fact  that  e~^  '^  /  0  for  all  oj  ,  it  is  clear  that  T  is 
one-to-one  on  any  space  of  Fourier  transformable  functions.  Further, 
since  the  inverse  of  the  multiplier  e*^  '^  has  no  Inverse  Fourier 
transform,  the  inverse  of  T  is  not  representable  as  a  convolution,  nor 
can  be  applied  to  the  general  space  of  all  Fourier  transformable 
functions.  Instead,  we  can  restrict  the  domain  of  T,  and  then 
respesent  its  inverse  as  a  convolution  on  the  range  of  T.  Many  such 
restricted  domains  are  possible.  In  the  next  section,  we  consider  T 
restricted  to  the  class  of  polynomials  of  degree  N  or  less. 
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III.   Polynomial  Domains 


Let  Pm  denote  the  space  of  polynomials  over  R  of  degree  less  than 
or  equal  to  N.  The  monomials  { 1 , x.x^ , . . . ,x^}  form  a  basis  for  V^  .  If 
this  basis  is  orthonormalized  with  respect  to  the  inner  product 


2 
(f,g)  =  J  f(x)  g(x)  e  '^   dx  ,  (3.1) 


then  the  basis   of   Hermite   polynomials   { H^^.H^ . . .  ,Hf^}   result.    The 
Hermites  can  be  represented  explicitly; 


V-)  = "'     ^  (-1)"  ?;„\,p' .  ^3-2^ 

n  -        m! (n-/m) ! 

m=0 


or  by  the  Rodrigues  formula; 


H  fx)  =  (-1)"  e^^  il  (e-'^S  ,  (3.3) 

"  dx" 


see,  e.g.,  [Courant  &  Hilbert]  or  [Lebedev]. 

Observation  1:  T  is  closed  on  P^  . 

Proof;  We  will  show  that  TH^  e  P^^  for  n  <  N. 


/^  (TH„)(y)  =  i  e-^y-""^      H„(x)  dx 


J  e-y^  e^^y  (-1)"  —  (e"''")  d 


dx" 


2y  I  e-y^  e^^^  (-1)""^  -^  (e'^^S  d 
■^  j„n-l 


•3-,<7 
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=  ^7  .  2y(THj^_i)(y) 


Thus 


(TH^)(x)  =  l^x""    .       •  (3. A) 


As   a   result   of   Observation  1 ,  T  is  an  isomorphisra  of  ?^  .   The 
inverse  of  T  on  P^.  is  clearly  given  by 

i=0        i=0 

Our  main  result  is  that  T"^  restricted  to  Pj^.  can  be   represented  by   a 
convolution  with  an  explicit  kernel  K^.(x): 

Theorera:  For  f  c  Pj^  and  g  =  Tf,  then 


f  =  S^  *  g 


(3.5) 


where 


2    ^^2      .k 
Kj,(x)  ^  e-'^      I      _ilii-- H2j,(x)  .  (3.6) 

k-O  ^   k!  2^ 

We  will  give  a  proof  below  using  direct  integration  (as  opposed  to 
using  Fourier  transform  distributions).  Note,  however,  that  K^,(x)  is 
not  the  unique  function  representing  T"^  on  P^  .  In  general,  the 
kernel  can  be  translated  by  any  function  which  yields  a  zero 
convolution  against  P^,  .    This   includes   all   functions  of  the  form 
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e   Hj^(x),  n  >  N.  The  stated  kernel  (3.6)  Is  unique  among  the  class   of 

2 
functions  of  the  form  e  ^   P(x),  where  P(x)  is  a  polynomial  of  degree  N. 

It   is   interesting   to   compare   the   form  of  Kj,(x)  with  standard 

enhancement  filters.   For  example,  for  N  =  3, 


K.(x)   =  i_  e"'^   (1  -  x2) 


nr 


2  dx2  /7 


1  „-x2  _  1  d^  (.Le-'^") 


Thus 


K3  *  g  =  [  1  -  J  -^-J  Tg  , 
dx*^ 


which  is  a  not  uncommon  high  emphasis  filter  (see,  e.g.,  the  papers  by 
E.  Machln  [Ratliff],  and  [Rosenfeld  and  Kak].  In  Figure  1,  we  display 
plots  of  Kvj  for  several  values  of  N. 

The  proof  of  the  theorem  depends  on  several  lemmas. 


Lemma  1 ;  ^n   ~     J 


i.  e"'^^  x"  dx  =  . 


,  n  odd 


/? 


nl 


2"(n/2)! 


n  even 


Proof:  For  n  odd,  e~^  x"  Is  an  Integrable  odd  function,  and  so  clearly 


Aj^  =  0.  For  n  =  2p,  p  >  1  , 


/7  A.   -  i  e"^  x2P  dj 


^2p 


-  i/  (-2x)  e-'^  x2p-l  d: 
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=  iPlL     J     e-'    x2p-2    dx   =  /7  iP;i 


'2p-2 


Since   Aq   =    1  , 


2p 


(2p-l)(2p-3)    ...    1      _      r2p) 


>m 


22"'p! 


The  formula  holds  for  p  >  0  . 


Lemma  2 : 


^2k,2p  =  I    ^^'"^     "Zk^'^)  ^^'^  dx  = 


(2p)l 


,   P  <  k 


22p-2k  (p_i,)! 


,   P  >  k 


Proof;  For  k  >  1 ,  p  >  1 , 


/V 


-2k,  2p 


y  e-^^  [(-1)21^  e-^i!i(e-^S]x2P  dx 
■^  dx^^ 


,2k 


=  i (e"'^  )  x2P  d: 


dx 


2k 


,2k-l 


-  2p  i  i (e-^  )  x2p-ldx 


dx 


2k- 1 


=  /V   (.2p)(2p-l)  C2j,_2,2p-2 


Clearly,  C2i^,^q  -  0,  k  >  1.  Using  Lemma  1,  Co^2p  '  (2p) ! /(Z^P/p ! ) ,  for  p 
>  0.  Combining,  C2j^^  2p  "  ^  ^°^   P  <  k,  and  for  p  >  k. 
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=2k,2p  =  (2p)(2p-l)  ...  (2p-2k+l)  .  Co^2p-2k 


(2p)!     (2p-2k)!    ^     (2p)! 
(2p-2k)!  22p-2k(p_^,),  '   22p-2k(p_k), 


Lemma  3 :  For  N  >  k, 


0  ,   k  odd 

'^^^    2'^(k/2)! 

Proof :  For  k  odd,   we   observe   from   (3.6)   that   Kv,(x)x    Is   an   odd 
integrable  function,  and  so  Integrates  to  zero.   For  k  =  2p, 


i  K^(x)x^dx  =  i  e"'^      Z    J"^^ H2i(x)  x^P  dx 

-«>  -«         1=0  '''tt  i !  2 


=   E 


^  (-l)i     (2p)I 


i=0i!2i  22p-2i(p-i)! 
=  ^2p)t  !*     P'    (_ni  d)P-i 


^'    (_  1)P  =  (-l)P   ^' 


2Ppl     2  2*^pl 

Proof-  of  the  Theorem;   By   equation   (3.4),   it   suffices   to  show  that 
Kj^  *  (2"x")  =  Hj^(x),  n  <  N.  We  have 

(K^  *  2"x")(y)   -  J  2"K^(x)(y-x)"dx 


00  n 

i  2%(x)  Z      (-1)^  [  {J  J  y""^  ^^   dx 
—         k=0 
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n 


=   I     .f"'  ,  (-1)^  •  y"-!^  .  d^^ 
k!  (n-k)  !  *^ 


k=0 


n/2    (_j)2m2n         (,^),    _2^ 

l!       i-   -rr: — r-n r — f:-  V~W 


m- 


=  0 


(2m)!(n-2m)!        yln^, 


n!      Z   .  .  "^^  ^,  2"-2"'  y"-2™  =  H^(y) 
-  m! (n-2m) ! 
m=0 


• 


The  theorem  above  could  have  been  proved  using  the  convolution 
theorem  and  by  computing  the  Fourier  transform  of  K^{x) .  We  will 
nonetheless  compute  V.^  in  order  to  show  that  the  multiplier  for  Kj^ 
approaches,  pointwise,  the  inverse  of  the  multiplier  for  the  operator  T 
(see  (2.2)). 

Observation  2:  Kj^(u) )  ■*■  e"^  '^   pointwise  as  N  -»•  » . 

N/2   .   .k       2 
Proof:  K„(co)  =  E   J-lii—.  rf  e"^  H2,^(x)]  (u) ,   where   F   stands   for   the 

k=0  -^  k!2'^  ' 
Fourier  tansform  operator.   Now, 


F[e-A2i,(x)J(a3)  =  f[  (-1)2^^  i^  (e"^  )]  M 


=   (io))2k/;r.  e^^/''  . 


Thus 


%'       (-^^"   .  (-1)^  .2^  /7  e-'/^ 
k=0  ♦^  1^1  2^ 


e 


-<ij 


9      N/2,2 
k-0 


^^2/4  ^^2/2  _  ^o)2/A 
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As   a   consequence   of   Observation  2,   we  see  that  K„(X)  does  not 
converge  polntwlse  to  any  function   as   N  ->   ^  ,   since   otherwise   the 

Fourier  transforni  of  that  function  would  be  e'^  '  ,  which  is  impossible. 

2   _  2 
Kj^(x)  does  converge  in  L  (e  '^  dx) ,  but  that  does   not   imply   polntwlse 

convergence   to   any   function.    We  accordingly  have  stable  deblurring 

when  using  the  kernels  Kj^(x),  where  stability  is  measured  in   terms   of 

9-2 
deviation   from   a   polynomial   of  degree  N,  and  the  L  (e  ^   dx)  norm  is 

used  as  the  metric. 
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IV.   Higher  Dimensions 

The  gaussian  blur  operator  is  given  by 


2 
Tf(x)  -  /  it""/2  e-(^-y)   f(y)  dy  .  (A.l) 


Due  to  the  separability  of  the  kernel  and  Fubini's  theorem,   T   can   be 
decomposed  into  n  iterated  blurrings: 


T  =  Tj  o  Tj  o  ...  o  T^  (A. 2) 

(T^f)(x)  =  ]  ±e~^''^~-^^'^    f^'^l  •••'  ^i  '  ••••  '^n)'^yi  (^-5) 

-oo  /tT 

Consider  a  polynomial  in  R*^: 


f(x)   =    I    a„x°  (^.^) 

|a  |<N 

a  =  ia^,...,a^)    ,  a^  e  Z  ,  a^  >  0  , 

a  1    a^ 
|a|=Iaj^,   X  =X]^...x^ 

For  fixed  x,  the  function  of  one  real  variable 

g(yj)   =   fCxj,...,  yj  ,  ....  Xj^) 

is  a  polynomial  of  degree  no  greater  than  N,  so 

Kj,  *  (Tg)  -  g  , 


'r -aI 
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where   T   is   the   standard   one   dimensional   blurring  operator  (2.1). 
Combining,  we  find  that 


f(x)  =  /   %(yi)  KN(y2)  •••  %(yn)  (Tf)(x-y)  dy        (A. 5) 


for  any  polynomial  f(x)  of  form  (A. 4).  Thus  deblurring  of  blurred 
polynomials  of  degree  N  can  be  accomplished  by  convolution  against  the 
kernel 

kJ(x)  =  Kf^Cxj)  •  K^(X2)  ...  K^Cx^).  (^-6) 

Thus   the   situation   in   higher   dimensions   is   similar   to   the   one 

dimensional   case.   The  deblurring  convolution  kernel  is  separable,  and 

2 
will  be  of  the  form  e~^  P(x),  where  P(x)  is  a  polynomial  of   degree   nN 

in  X  e  R".   Figure  2  shows  a  plot  of  Kj^(x)  for  n  =  2,  N  =  3. 
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V.  Some  Experimental  Results  and  Comments 

In  Figure  3,  we  show  a  512  by  512  pixel  array  digitized  to  256 
gray  levels,  focused  and  imaged  as  carefully  as  possible.  The  two 
dimensional  blurring  operator  T  was  applied  to  obtain  Figure  A,  using 
an  Interplxel  distance  of  0.2  in  both  dimensions.  A  VICOM  Systems, 
Inc.,  image  processing  computer  was  used  to  perform  the  convolution, 
using  12  bits  of  significance  in  the  calculation.  The  convolution  was 
performed   in  each  dimension  separately,  since  the  kernel  is  separable, 

using   coefficients   obtained   by   integrating   0.2   intervals   of   the 

2 
gaussian  function  1 //JT  (e"'^  ).   In  a  similar  manner,  deblurring  kernels 

kJJ^Cx)  were  applied   to   obtain   Figures   5,   6   and   7.   Again,   a   0.2 

interplxel  distance  was  used,  so  that  the  deblurring  extent  matched  the 

blurring  operator. 

We  observe  that  as  N  Increases,  successively  better  deblurrings 
are  obtained.  In  no  case,  however,  is  the  reconstruction  visually 
Indistinguishable  from  the  original  (Figure  3).  Further,  as  N 
increases,  the  fluctuations  in  K„(x)  become  more  violent,  and  will  lead 
to  numerical  inaccuracies  in  the  deblurring  process. 

This  inabiity  to  accurately  reconstruct,  and  the  artifacts  that 
appear  as  a  result  of  the  deblurring  process,  reflect  the  fact  that 
natural  Image  data  Is  not  represented  well  locally  by  Nth  order 
polynomials.  This  does  not  preclude  the  possibility  that  some 
Invertible  transformation  of  the  image  data  j^  well  represented  locally 
by  polynomials,  and  so  will  perform  better.  Classes  of  signals  other 
than  natural  Images  may  admit  better  approximations  by  polynomials. 
Perhaps   a   treatment   of   the  same  deblurring  problem  using  a  class  of 
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analytlc  functions  different  than  Nth  order   polynomials   Is   possible, 
and  will  lead  to  better  performance  on  Image  data. 

The  real  issue  of  this  paper  Is  the  search  for  an  Intermediate 
representation  of  signal  data  In  a  fashion  that  preserves  only  that 
information  essential  to  the  interpretation  of  the  data.  For  imagery, 
this  means  that  reconstruction  of  a  visually  similar  picture  should  be 
possible  from  a  good  representation.  This  paper  has  shown  that  a 
blurred  version  of  data  is  a  suitable  representation  when  the  original 
data  is  well  approximated  locally  by  polynomials  of  bounded  degree. 
While  this  may  not  describe  general  Images,  it  may  nonetheless  be  the 
key  to  alternate  representations  involving  gaussian  blurrlngs. 


The  supervision  of  this  work  by  Professor  Steven  Zucker  of  McGill 
University  is  gratefully  acknowledged. 
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Figure  1.   One  dimensional 
deblurring  kernels  K^(x), 
N=3,5,9. 
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Figure  2.  Two  dimensional  deblurring  kernel,  N  -  3, 
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Figure  3:   Original  image,  the  Saltair  resort,  circa  1930. 
(Photo  from  the  Utah  State  Historical  Society.) 
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Figure  4:   Blurred  version  of  the  original  image  in 
Figure  3 . 


-23- 


Figure  5:   Deblurring  kernel  K  applied  to  the  blurred 
image  in  Figure  4 . 


-24- 


Figure  6:   Deblurring  kernel  K  applied  to  the  blurred 
image  in  Figure  4. 
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Figure  7:   Deblurring  kernel  K   applied  to  the  blurred 
image  in  Figure  4 . 
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